#!/bin/bash
set -e

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------

# java -jar $gatk \
# -R $ref_genome \
# -T SelectVariants \
# -V $1 \
# -o $2 \
# -select "DP > 30" \
# -select "AF > 0.2" \
# -selectType SNP


# java -jar $gatk \
# -R $ref_genome \
# -T SelectVariants \
# -V $1 \
# -o $2 \
# -select 'vc.getGenotype("ca50mor1").getAD().1  / vc.getGenotype("ca50mor1").getDP() >  0.21' \
# -selectType SNP

# -select "DP <= 30" \
# -select 'vc.getGenotype("NA12878").getAD().0 > 10' \


# java -Xmx$java_memory -jar $gatk \
	# -R $ref_genome \
	# -T VariantFiltration \
	# -o $2 \
	# --variant $1 \
	# --filterExpression "FS > 60.0" \
	# --filterName "FS_more_than_60" \
	# --filterExpression "MQ < 40.0" \
	# --filterName "MQ_less_than_40" \
	# --filterExpression "HaplotypeScore > 200.0" \
	# --filterName "HaplotypeScore_more_than_200" \
	# --filterExpression "af < 0.2" \
	# --filterName "Allele_frequency_less_than_20%" \
	# --filterExpression "DP < 30" \
	# --filterName "Depth_less_than_30"
	
echo
echo
echo VEP annotation no pick
variant_effect_predictor.pl \
--offline \
--force_overwrite \
--dir $data_path/ensembl/ \
--fasta $ref_genome \
--refseq \
--hgvs \
--format vcf \
--vcf \
--input_file 1.dbnsfp.snpeff2gatk.vcf \
--output_file 1.vep.vcf


cat 1.vep.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $1","$2","$3","$4}else{print $0}}' \
> 1.header_modified.vcf



echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift filter 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
filter \
"(FILTER = 'PASS')" \
1.header_modified.vcf \
> 1.high.vcf \
2>>$log


# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo $snpsift filter 2>>$log 1>&2
# time \
# java -Xmx$java_memory -jar $snpsift \
# filter \
# "(FILTER != 'PASS')" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 2>>$log



$bgzip -c 1.header_modified.vcf > 1.header_modified.vcf.gz
$tabix 1.header_modified.vcf.gz

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo bcftools view for filter FILTER 2>>$log 1>&2
time \
$bcftools view \
-f "FP_VARIANTCALLER,FP_PRIMERDIMER" \
1.header_modified.vcf.gz \
| bgzip -c \
> 1.fp.vcf.gz \
2>>$log

gunzip -c 1.fp.vcf.gz > 1.fp.vcf

$tabix 1.fp.vcf.gz

$bgzip -c 1.high.vcf > 1.high.vcf.gz
$tabix 1.high.vcf.gz
$bcftools concat -aD 1.high.vcf.gz 1.fp.vcf.gz |bgzip -c > 1.high_fp.vcf.gz
$tabix 1.high_fp.vcf.gz

if [ ! -d isec ]; then
	mkdir isec
fi
$bcftools isec -p ./isec 1.header_modified.vcf.gz 1.high_fp.vcf.gz

cp ./isec/0000.vcf 1.low.vcf